rm(list=ls())


# Check that required packages are installed:
want = c("foreign", "ggplot2", "MASS", "reshape2", "plyr", "dplyr", "ggrepel", 
         "lme4", "coda", "relaimpo", "relimp", "arm", "boot", "entropy", "here")
have = want %in% rownames(installed.packages())
if ( any(!have) ) { install.packages( want[!have] ) }
# load packages
junk <- lapply(want, library, character.only = TRUE)
rm(have,want,junk)

options(scipen = 999)

#=============================================================================================================
# Load Image
load("allmodels.RData")
#=============================================================================================================


#=============================================================================================================
# RELATIVE IMPORTANCE
#=============================================================================================================
#---------------#
# Grömping 2007 #
#---------------#
n.sst <- c("clasd1","clasd2","clasd4","clasd5","standard_gsd","standard_gsd_sq",
           "churchat_gsd","churchat_gsd_sq","relig_gsd","relig_gsd_sq")
n.iss <- c("i2ImmCus_gsd","i1PriEnt_gsd","i3SamSex_gsd","i1StaOwn_gsd","i3AboFree_gsd","i1IntEco_gsd",
           "i3HarSen_gsd","i1OrdRed_gsd","i3ObeAut_gsd","i3WorCut_gsd","i2ImmDec_gsd","i2ImmCus_gsd_sq",
           "i1PriEnt_gsd_sq","i3SamSex_gsd_sq","i1StaOwn_gsd_sq","i3AboFree_gsd_sq","i1IntEco_gsd_sq",
           "i3HarSen_gsd_sq","i1OrdRed_gsd_sq","i3ObeAut_gsd_sq","i3WorCut_gsd_sq","i2ImmDec_gsd_sq")
rsq <- cbind(unique(data.un$syslab),
             data.frame(matrix(NA, nrow = length(unique(data.un$syslab)), ncol = 3)))
names(rsq) <- c("syslab","sst","iss","part")
for (i in unique(data.un$syslab)) {
  pties <- names(table(subset(data.un, syslab == i)$close_party_full))
  ivs <- c("female", "age", "edlow", "edhigh", n.sst, n.iss)
  np <- paste0("close_str_", pties)
  fml <- as.formula(paste0("lrsp ~ ", paste(c(ivs, np), collapse = "+")))
  no.na.data <- subset(data.un, syslab == i)
  no.na.data <- no.na.data[rowSums(is.na(no.na.data[, np])) == 0, ]
  rsq[rsq$syslab == i, 2:4] <- calc.relimp(formula = fml,
                                           data = no.na.data,
                                           groups = list(n.sst, n.iss, np))$lmg[c("G1", "G2", "G3")]
  rm(no.na.data)
  print(i)
}

#=============================================================================================================
# Plots
#=============================================================================================================
# Plot of all R-Squares
rsq.melt <- melt(rsq, id.vars = c("syslab"))
rsq.melt$IS <- ifelse(rsq.melt$variable == "sst", "Social Structure",
                      ifelse(rsq.melt$variable == "iss", "Issues", "Partisanship"))
rsq.melt$IS <- factor(rsq.melt$IS, levels = c("Issues", "Social Structure", "Partisanship"))
rsq.melt <- data.frame(rsq.melt %>% group_by(syslab) %>% dplyr::mutate(rsq.mean = mean(value)))
rsq.melt$syslab <- factor(rsq.melt$syslab, levels = unique(rsq.melt[order(rsq.melt$rsq.mean), "syslab"]))


quartz(type = 'pdf', file = 'rsquares_gromping.pdf',
       width = 10, height = 6, dpi=300)
ggplot(rsq.melt, aes(y = value, x = syslab, group = IS, fill = IS)) + 
  geom_bar(position = "dodge", stat = "identity") + 
  scale_fill_manual(values = c("gray70", "gray50", "gray30"),"") +
  theme_bw() + 
  theme(legend.position="top") +
  xlab("") + ylab("R-Squares")
dev.off()


# Compute ratio
rsq$ratio <- log(rsq$part / rowSums(rsq[, c("iss", "sst")]), 2)
# Plot ratio
quartz(type = 'pdf', file = 'ratio_gromping.pdf',
       width=7,height=5, dpi=300)
ggplot(rsq, aes(x = syslab, y = ratio)) + 
  geom_point() + 
  coord_flip() +
  scale_x_discrete(limits = (rsq[rev(order(rsq$ratio)), "syslab"])) +
  theme_bw() +
  xlab("") + 
  ylab("Importance of Party Attachment over Substantive\nConsiderations for Left-Right self-placement (ratio)")
dev.off()

# Compare with ratio obtained with Silber et al. method
ratio.comp <- merge(rsq[, c("syslab", "ratio")], ratio.lm.boot[, c("syslab", "ri_l2_m")])


quartz(type = 'pdf', file = 'methods_comparison.pdf',
       width = 7, height = 6, dpi=300)
ggplot(ratio.comp, aes(x = ri_l2_m, y = ratio)) +
  #geom_abline(aes(intercept = 0, slope = 1), col = "gray50") +
  geom_smooth(method = "lm", se = F) +
  geom_point() + geom_text_repel(aes(label = syslab)) +
  xlab("Silber et al. (1995) method") +
  ylab("Grömping (2007) method") +
  ggtitle(paste0("Correlation = ", round(with(ratio.comp, cor(ri_l2_m, ratio)), 3))) +
  scale_x_continuous(limits = c(-1.22, 2.26)) +
  scale_y_continuous(limits = c(-1.22, 2.26)) +
  theme_bw()
dev.off()


